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Abstract 

We investigate the discretization of Darcy flow through fractured porous media on 
general meshes. We consider a hybrid dimensional model, invoking a complex network of 
planar fractures. The model accounts for matrix-fracture interactions and fractures acting 
either as drains or as barriers, i.e. we have to deal with pressure discontinuities at matrix- 
fracture interfaces. The numerical analysis is performed in the general framework of 
gradient discretizations which is extended to the model under consideration. Two families 
of schemes namely the Vertex Approximate Gradient scheme (VAG) and the Hybrid Finite 
Volume scheme (HFV) are detailed and shown to satisfy the gradient scheme framework, 
which yields, in particular, convergence. Numerical tests confirm the theoretical results. 
Gradient Discretization; Darcy Flow, Discrete Fracture Networks, Finite Volume 


1 Introduction 

This work deals with the discretization of Darcy flows in fractured porous media for which 
the fractures are modelized as interfaces of codimension one. In this framework, the d — 1 
dimensional flow in the fractures is coupled with the d dimensional flow in the matrix leading 
to the so called, hybrid dimensional Darcy flow model. We consider the case for which the 
pressure can be discontinuous at the matrix fracture interfaces in order to account for fractures 
acting either as drains or as barriers as described in [10], [12] and [3]. In this paper, we will 
study the family of models described in [12] and [3]. 

It is also assumed in the following that the pressure is continuous at the fracture intersec¬ 
tions. This corresponds to a ratio between the permeability at the fracture intersection and 
the width of the fracture assumed to be high compared with the ratio between the tangential 
permeability of each fracture and its length. We refer to [14] for a more general reduced model 
taking into account discontinuous pressures at fracture intersections in dimension d = 2. 

The discretization of such hybrid dimensional Darcy flow model has been the object of 
several works. In [10], [11], [3] a cell-centered Finite Volume scheme using a Two Point Flux 
Approximation (TPFA) is proposed assuming the orthogonality of the mesh and isotropic per¬ 
meability helds. Cell-centered Finite Volume schemes have been extended to general meshes 
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and anisotropic permeability fields using MultiPoint Flux Approximations (MPFA) in [13], [16], 
and [2]. In [12], a Mixed Finite Element (MFE) method is proposed and a MEE discretiza¬ 
tion adapted to non-matching fracture and matrix meshes is studied in [6]. More recently the 
Hybrid Einite Volume (HFV) scheme, introduced in [8], has been extended in [27] for the non 
matching discretization of two reduced fault models. Also a Mimetic Einite Difference (MFD) 
scheme is used in [1] in the matrix domain coupled with a TPFA scheme in the fracture network. 
Discretizations of the related reduced model [28] assuming a continuous pressure at the matrix 
fracture interfaces have been proposed in [28] using a MFE method, in [20] using a Control 
Volume Einite Element method (CVEE), in [19] using the HEV scheme, and in [19, 5] using an 
extension of the Vertex Approximate Gradient (VAG) scheme introduced in [7]. 

In terms of convergence analysis, the case of continuous pressure models at the matrix frac¬ 
ture interfaces [28] is studied in [19] for a general fracture network but the current state of 
the art for the discontinuous pressure models at the matrix fracture interfaces is still limited 
to rather simple geometries. Let us recall that the family of models introduced in [12] and [3] 
depends on a quadrature parameter denoted by G [|, 1] for the approximate integration in 
the width of the fractures. Existing convergence analysis for such models cover the case of one 
non immersed fracture separating the domain into two subdomains using a MFE discretization 
in [12] or a non matching MEE discretization in [6] for the range ^ G (^,1]- [3], the case of 

one fully immersed fracture in dimension d = 2 using a TPEA discretization is analysed for the 
full range of parameters ^ G [|, !]• 

The main goal of this paper is to study the discretizations of such models and their con¬ 
vergence properties by extension of the gradient scheme framework. The gradient scheme 
framework has been introduced in [7], [22], [21] to analyse the convergence of numerical meth¬ 
ods for linear and nonlinear second order diffusion problems. As shown in [22], this framework 
accounts for various conforming and non conforming discretizations such as Finite Element 
methods. Mixed and Mixed Hybrid Finite Element methods, and some Einite Volume schemes 
like symmetric MPEA, the VAG schemes [7], and the HFV schemes [8]. 

Our extension of the gradient scheme framework to the hybrid dimensional Darcy flow model 
will account for general fracture networks including fully, partially and non immersed fractures 
as well as fracture intersections in a 3D surrounding matrix domain. Each individual fracture 
will be assumed to be planar. The framework will cover the range of parameters ^ G (^,1] 
excluding the value ^ in order to allow for a primal variational formulation. 

Two examples of gradient discretizations will be provided, namely the extension of the VAG 
and HFV schemes defined in [7] and [8] to the family of hybrid dimensional Darcy flow models. 
In both cases, it is assumed that the fracture network is conforming to the mesh in the sense 
that it is defined as a collection of faces of the mesh. The mesh is assumed to be polyhedral 
with possibly non planar faces for the VAG scheme and planar faces for the HFV scheme. Two 
versions of the VAG scheme will be studied, the first corresponding to the conforming Pi finite 
element on a tetrahedral submesh, and the second to a finite volume scheme using lumping for 
the source terms as well as for the matrix fracture fluxes. The VAG scheme has the advan¬ 
tage to lead to a sparse discretization on tetrahedral or mainly tetrahedral meshes. It will be 
compared to the HFV discretization using face and fracture edge unknowns in addition to the 
cell unknowns. Note that the HFV scheme of [8] has been generalized in [23] as the family of 
Hybrid Mimetic Mixed methods which which encompasses the family of MFD schemes [24]. In 
this article, we will focus without restriction on the particular case presented in [8] for the sake 
of simplicity. 

In section 2 we introduce the geometry of the matrix and fracture domains and present 
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the strong and weak formulation of the model. Section 3 is devoted to the introduction of the 
general framework of gradient discretizations and the derivation of the error estimate 3.3. In 
section 4 we dehne and investigate the families of VAG and HFV discretizations. Having in 
mind applications to multi-phase flow, we also present a Finite Volume formulation involving 
conservative fluxes, which applies for both schemes. In section 5, the VAG and HFV schemes 
are compared in terms of accuracy and GPU efficiency for both Gartesian and tetrahedral 
meshes on hererogeneous isotropic and anisotropic media using a family of analytical solutions. 

2 Hybrid dimensional Darcy Flow Model in Fractured 
Porous Media 

2.1 Geometry and Function Spaces 

Let G denote a bounded domain of IR'^, d = 2,3 assumed to be polyhedral for d = 3 and polyg¬ 
onal for d = 2. To hx ideas the dimension will be hxed to d = 3 when it needs to be specihed, 
for instance in the naming of the geometrical objects or for the space discretization in the next 
section. The adaptations to the case d = 2 are straightforward. 

Let F = UjgjFj and its interior F = F \ 9F denote the network of fractures F* C G, f G /, 
such that each Fj is a planar polygonal simply connected open domain included in a plane Vi 
of It is assumed that the angles of Fj are strictly smaller than 27 r, and that Fj fl F^ = 0 for 
all i ^ j ■ 

For all i & I, let us set Sj = SFj, with ns. as unit vector in Vi, normal to Sj and outward 
to Fj. Further Sjj = Ej fl S^, j e I \ {i}, Sj,o = Sj fl 5G, Sj^jv = Sj \ (Ujg/\p} U Sj^o), 

E = [j{i,j)(^ixi,i^ji^i,j \ ^i,o) and Eq = Uie/ ^*.0- is assumed that Ej^o = Tj fl dn. 




Figure 1: Example of a 2D domain D and 3 intersecting fractures Fj, z = 1,2, 3. We might dehne 
the fracture plane orientations by = ai,a“(l) = for Fi, q;+( 2) = ai,Q!“(2) = 0^2 for 

F2, and Q!’''(3) = a3,Q!“(3) = 0^2 for F3. 


We will denote by dr(x) the d — 1 dimensional Lebesgue measure on F. On the fracture 
network F, we dehne the function space L^(F) = {n = {vi)i^i,Vi G L^(Fj),z G /}, endowed 
with the norm ||n||L2(r) = (Xlie/llD|li2('r-))^ and its subspace id^(F) consisting of functions 
V = {vi)i^i such that Vi G dd^(Fj), z G / with continuous traces at the fracture intersections 
Ejj, j E I \ {z}. The space H^iT) is endowed with the norm ||n||zzi(r) = {JZiei llull^pr.))^- 
We also dehne it’s subspace with vanishing traces on Eq, which we denote by 

On D\F, the gradient operator from id^(f2\F) to L‘^{QY is denoted by V. On the fracture 
network F, the tangential gradient, acting from id^(F) to L^(F)'’*“^, is denoted by V,-, and such 
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that 


Vt-W 

where, for each i E I, the tangential gradient Vn is dehned from to by hxing 

a reference Cartesian coordinate system of the plane Vi containing Fj. We also denote by div^-. 
the divergence operator from H^^y(Ti) to L‘^(Ti). 

We assume that there exists a hnite family (rQ,)aex all a G y holds: Fq, C F 

and there exists a lipschitz domain C hl\F, such that F^ = (9a;Q,nF. For a E x and an apro- 
priate choice of C / we assume that Fa = U*e/c r*- Furthermore should hold F = F„. 

We also assume that each Fj C F is contained in Fq, for exactly two a E x and that we can 
dehne a unique mapping i \ — > (Q;+(i), a~{i)) from / to y x y, such that F* C F^+p) flFo-p) and 
7^ a~{i) (cf. hgure 1). For all i E I, Q;^(i) dehnes the two sides of the fracture Fj in r2\F 
and we can introduce the corresponding unit normal vectors na±{i) at Fj outward to Ci;Q,±(j), such 
that no+(j) + = 0. We therefore obtain for a G y and a.e. x G Fq, a unique unit normal 

vector nQ(x) outward to cUq. A simple choice of (FQ)Qg^ is given by both sides of each fracture 
i E I but more general choices are also possible such as for example the one exhibited in hgure 1. 

Then, for a G y, we can dehne the trace operator on Fq: 

7Q://'(ff\r)^L2(FQ), 
and the normal trace operator on Fq outward to the side a: 

7„,Q:i/div(^\n^^^'(r«). 

We now dehne the hybrid dimensional function spaces that will be used as variational spaces 
for the Darcy how model in the next subsection: 

V = ff\Q\T) X iF^(F), 

and its subspace 

D° = iF^^(D\r)xFf^^(F), 

where (with jgQ: iF^(D\F) —)■ L‘^{dQ) denoting the trace operator on dQ) 

\ F) = {n G H\n\T) I Xdnv = 0 on 9D}, 

as well as 

W = WmX Wf, 

where 

Wm = {qm e hfdiv(^ \ I 7n,aqm ^ T^(Fq) for all a G y} and 
^/ = {q/ = (q/,*)*e/ I q/,i e iF^jy(Fi) for alH G / 



On V, we dehne the positive semidehnite, symmetric bilinear form 

{{u^,Uf),{v^,Vf))v = / Vu^-Vvmd^+ / VrUfVrVfdri^) 

Jn Jr 

+ ^/ {laUm-Uf){'yaVni-Vf)dT{x.) 

aex 
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for {um,Uf), {vm,Vf) G V, which induces the seminorm \{vm,Vf)\v- Note that (•, ■)v is a scalar 
product and | ■ |y is a norm on V^, denoted by || ■ Hyo in the following. 


We dehne for all (pm,P/), (Qm^Q/) ^ W the scalar product 

((Pm,P/), (qm,q/))iy = [ Pmqmdx+ [ divp^ ■ divq^dx 


y'p/q/dr(x) + J div^pf ■ div^qf dr{:} 

/ (7n,aPm ■ 7n,aqm)dr(x), 

^,^-vx ^rrv 


which induces the norm ||(qm, q/)||w; and where we have used the notation div^-pj = div^-.p/^i 
on Fj for alH G / and p/ = (p/,i)jg 7 G Wf. 

Using similar arguments as in the proof of [15], example II.3.4, one can prove the following 
Poincare type inequality. 


Proposition 2.1 The norm || ■ ||yo satisfies the following inequality 

ll'^m|lHi(0\r) + llwll^^hr) < Cp\\{Vm,Vf)\\vO, 

for all {vm,Vf) G U°. 


( 1 ) 


Proof We apply the ideas of the proof of [15], example II.3.4 and assume that the statement 
of the proposition is not true. Then we can dehne a sequence {vi)i^u in U°, such that 


\vi\\vo < y. 


( 2 ) 


||nz||i^i = 1 and 
where, for this proof, || ■ ||//i = || ■ + II ' ll^qr)- The imbedding 

II ■ IIhO X T^(r), II ■ ||l2(q) + II ■ ||L2(p)j 

is compact, provided that r2\r has the cone property (see [18], theorem 6.2). Thus, there is a 
subsequence of (nz)zeiisi and v G L‘^{Q) x T^(r), such that 

Vfj^ —)■ V in L^(r2) X L^{T). 

On the other hand it follows from (2) that 

—)■ 0 in L^(0) 

Since (U°, || ■ ||jpi) is complete, we have 


with 


in U°, 


|n||yo = lim ||n^||yo = 0. 


Since || ■ ||yo is a norm on U°, we have n = 0 G U°, but ||n|| = 1, which is a contradiction. □ 
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Remark 2.1 With the precedent proof it is readily seen that inequality (1) holds for all func¬ 
tions V E V whose trace vanishes on a subset of d{Q\T) with positive surface measure. The 
requirement is that v has to be in a closed subspace of (V, || ■ \\h^) for which || ■ H^o is a well 
defined norm. 

The convergence analysis presented in section 4 requires some results on the density of 
smooth subspaces of V and W, which we state below. 

Definition 2.1 1. Cq is defined as the subspace of functions in \ T) vanishing on a 

neighbourhood of the boundary dTl, where \ T) C \ T) is the set of functions 

ip, such that for all is. E Q there exists r > 0, such that for all connected components u of 
{x + y G [R*^ I |y| < r} n (12 \ T) one has p E C°°{u). 

2. = 7 r(C“( 12 )) is defined as the image o/C“(12) of the trace operator'yr- H^i^) —t 
L^T). 

3. C^^ = C^{Q\Tf. 

4- = {q/ = (q/,i)ie/ I q/,i e = o on S, q/,* ■ ns, = 

0 on Sj^AT, i E I}. 

Let us hrst state the following Lemma that will be used to prove the density of x in 
W. 

Lemma 2.1 LetVm G 2^^(12), Vf E L^(r), G E L‘^{ytY, H E L‘^{VY~^ and Ja E //^(Lq), a E x 
such that 

/ {G-(:\rn + nmdivq„^)dx + / (iL ■ q/ + n/div^q/)dr(x) + V] / 7n,aqmdr(x)(E - O/) = 0 

Jr dr. 

(3) 


for all (qm,q/) e G^^ x G^^. Then holds {vm,Vf) E V^, {G,H) = {Vvm,VrVf) and = 
Vf - XaVm for a Ex- 

Proof Firstly, for all q^ E C“(12\r)‘^, we have 

/ (G ■ q^ + nmdivq„^)dx = 0 
Jn 

and therefore Vm E H^{Q\T) and Vvm = G. 

For a.e. x G (912, there exists an open planar domain u CC (912\(9F containing x such that 
for all / G G“(a;) there exists q^ E G^^ with 


bnanq m 


7n,aqm 


J / on to, 

I 0 on (912\a;, 

0 on Fo, a G X, 


where Xngn denotes the normal trace operator on the boundary of 12. From (3), taking q/ = 0, 
we obtain 


0 = / (Vom ■ q,„ + nmdivq^)dx = / XdnVm'Jnen^mdrix) = / XdnVmfdrix). 

Jo. Jon Juj 
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where 7 ^^ denotes the trace operator on the boundary of We deduce 'yan'Vm = 0 a.e. on 
dn\dT. Hence G 

Further, for a.e. x G Fq, there exists an open planar domain uja CC Fq containing x such 
that for all g G C“(a;Q,) there exists G with 


'yn.,a^m 
Fn,/3Qm 

Fna^q m 


( g on cuo, 

\ 0 on Fo\a;a, 

0 on F^, for (3 a, 
0 on (9H. 


From (3) we obtain 

0 = / (VUm ■ qm + nmdivqm)dx + ^ / 'yn,aqm{Ja - Vf)dT{x.) 

Jn 

= / “ W + 7a^'m)dr(x) = / g{Ja -Vf+ 7anm)dr(x). 

V r Q; V 

We deduce Ja = Vf — 'yaVm a.e. on Fq,, q; G y. 

Next, for all q/ G C“(Fj)‘^“^, i e I, we have from (3) 


/ {H ■ qf + n/divq/)dr(x) = 0 

Jr, 

and therefore u/l'p, G H^{Ti) for i G / and Vri^/fri = H\r,- 

Let i,j E I, i j. For a.e. x G Sjj \ Sj^o there exists an open interval CC Sjj \ Sj^o 

containing x such that for all h G C^(cij) there exists s G with 

h 7ns Cj) 

^ 3 

Ins^s = 0 on Sfc\cy, ke I. 

From (3) we obtain 


0= (V rVf ■ s + n/divT-s)dr(x) 



7sW/)7ns^sda(x), 


d(T(x) denoting the d — 2 dimensional Lebesgue measure on S. We deduce 'yj:,Vf = 7 e W a.e. 
on Tiij \ Sjo, i,j G I,i 7 ^ j. The proof of 7soW = 0 a.e. on Sq goes analogously. Hence 
■Of € Hun □ 


Proposition 2.2 Cq X is dense in V^. 


Proof Firstly, note that we have 

^11 VMm||E2(Q)<i + II VT-'lt/||j;^2(-r)ti-i^ < II (rim, M/) II V'O 

— ■ ^11 Vrtm 11^2(0)'* + II ^rW/ll , 

i.e. Il'llyo is equivalent to the standard norm || V ■ ||^ 2 (Q)d + || ■ ||j;^ 2 ('p)d-i on The density of 

Cq in Hq^{Q\T) being a classical result, we are concerned to prove the density of in H^^(T) 
in the following. Since C 7r(Lfd(H)), we can dehne = 'y^^{H^^{T)) C iLd(^)- 

Proposition 2 of [19] it is shown that C^{Q) is dense in (H°, ||V ■ ||l2(o)'* + l|Vr7r ■ ||L2(r)‘*-i)- 
Hence is dense in ||Vt- ■ ||p 2 ('r)d-i). □ 
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Proposition 2.3 x is dense in W. 

Proof Since IP/ is a closed subspace of the Hilbert space riie/ -^div(^O) linear form I G W'^ 
is the restriction to Wf of a linear form still denoted by I in Then, for some 

/ e L\r) and g G holds 


< ^, q/ >= X] / (^ ■ q/ + / ■ div^q/)dr(x), 

i&l 

for all q/ G Wf. Let us assume now that < /, (^ >= 0 for all G Corresponding to 

Lemma 2.1 holds / G From the dehnition of Wf we conclude that < I, q/ >= 0 for all 

q/ ^ 

Let now / G IFj^. Then there exist / G g G and ha G L‘^{Ta) {a G y), such 

that 

ha7n,aqmdr(x), 

for all qm G ITm. Furthermore, let us assume that <l,<~p >= 0 for all G From Lemma 

2.1 we deduce that / G Hq^{VL \ F), that g = V/ and that ha = Jaf (« ^ x)- Using this, we 

conclude, again by the rule of partial integration, that < I, q^ >= 0 for all q^ G IFm. □ 

2.2 Single Phase Darcy Flow Model 

2.2.1 Strong formulation 

In the matrix domain H\F, let us denote by G the permeability tensor such that 

there exist Am > > 0 with 

Amicr < (A„.(x)C,C) < Amici' for all C e R",x G H, 

Analogously, in the fracture network F, we denote by A/ G the tangential 

permeability tensor, and assume that there exist A/ > A/ > 0, such that holds 

A,Id" < (A/(x)C,C) < A;|C|" tor all C £ R^-dx e T. 

At the fracture network F, we introduce the orthonormal system 
(ti(x), T 2 (x), n(x)), dehned a.e. on F. Inside the fractures, the normal direction is assumed 
to be a permeability principal direction. The normal permeability A/,n £ L°°(F) is such that 
A/^n — A/,n(x) < A/,n for a.e. x G F with 0 < A/.n — A/,n- We also denote by df G L°°(T) the 
width of the fractures assumed to be such that there exist df > df > 0 with 


< /, qm >= / ig-qm + f ■ divq^ 


dx£^/ 

a&X ^ 


df < d/(x) < df 

for a.e. x G F. Let us dehne the weighted Lebesgue d — 1 dimensional measure on F by 
dTfipd) = (i/(x)dr(x). We consider the source terms hm G (resp. hf G L^(F)) in the 

matrix domain H \ F (resp. in the fracture network F). The half normal transmissibility in the 

2 A 

fracture network is denoted by T/ = 



Given ^ G (|, 1], the PDEs model writes: find {um,Uf) G (qm,q/) G W such that: 


div(qm) 

qm 


< 


diVr,(q/) - 7n,a+(i)qm 


7n,«±(i)qm 

7n,«“(i)qm 

q/ 


rp 

2^—1 i.^'y “1“ 

dfhf 

-df AfVrUf 


on Q\T, 
on G \ r, 

Uf) on hi, i G I, (4) 
on Fj, i ^ 
on r, 


2.2.2 Weak formulation 

The hybrid dimensional weak formulation amounts to find {um,Uf) G satisfying the following 
variational equality for all {vm,Vf) G V^: 


AmVUra ' Vu^dx + / AfVrUf ' VrVfdiTf{^) 


E 

i&I 


Tf 


T 

A I 


2^-1 


+ (1 - 07/3Wm - M/) (laVm “ V/) dr(x) (' 5 ) 


(a,/3)e{(o±(i),oT(i))} 


- / hmt'mdx- / h/t;/dr/(x) = 0. 

Jo Jr 

The following proposition states the well posedness of the variational formulation (5). 

Proposition 2.4 For all ^ G (|, 1], the variational problem (5) has a unique solution {um,Uf) G 
which satisfies the a priori estimate 


Ujm Uf)\\vo < C[ \\hm\\L^(n) + ||h/||L2(r)), 


with C depending only on Cp, \^, \f, df, df, and Ayn- addition (qm,q/) = 

— {AmVumidfAfVrUf) bclongs to W. 

Proof Using that for all G (|, 1] and for all (a, b) G one has 

+ (1 — 0 ^)® + + (1 “ 0 ®)^ — 2 ^ — 1 

the Lax-Milgram Theorem applies, which ensures the statement of the proposition. □ 


3 Gradient Discretization of the Hybrid Dimensional Model 

3.1 Gradient Scheme Framework 

A gradient discretization V of hybrid dimensional Darcy fiow models is defined by a vector space 
of degrees of freedom Xp = Xp^ x Xp ^, its subspace satisfying ad hoc homogeneous boundary 
conditions X^ = X^^ x Xl^^, and the following gradient and reconstruction operators: 

• Gradient operator on the matrix domain: Vp^ : Xp^ —)■ Lfi{VLY 

• Gradient operator on the fracture network: Xpj : Xp^ —)■ 

• A function reconstruction operator on the matrix domain: 

np^ : ^ 
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• Two function reconstruction operators on the fracture network: 

^ L2(r) and ^ ^^(r) 

• Reconstruction operators of the trace on Tq for a E X'- 

n?,„ : Ap„ ^ L2(r„). 

The space X® is endowed with the seminorm 

1 

aSx 

which is assumed to dehne a norm on X^. 

The following properties of gradient discretizations are crucial for the convergence analysis 
of the corresponding numerical schemes: 


Coercivity: Let P be a gradient discretization and 

\\^Vm'^v^\\L2{n) + W'n.VfVVfWi'^ir) 


Cd = max 


\\{vv^,VVf)\\v 


A sequence (T’*)zeiisi of gradient discretizations is said to be coercive, if there exists Cp > 0 such 
that Cpi <Cp for all I G N. 


Consistency: Let T* be a gradient discretization. Torn = {um,Uf) G C^anduD = {vp^,vp^) G 
X^ let us dehne 


s{vp,u) = \\Xv^Vv^ — V-Um||L 2 (Q)d + ||l, 2 (r)d-i 

+ ~ Um\\L2{n) + \\Tlvj>vp^ - -U/I|x 2 (r) 

+ WUp^w^ - n/||L2(r) + Zlaex “ XaUmWi^ir^)- 

and Spin) = min„^gxo s{vp,u). A sequence (T’^)«giH of gradient discretizations is said to be 
consistent, if for all u = {um,Uf) G holds 

lim Spi{u) = 0. 

1^00 


Limit Conformity: Let "D be a gradient discretization. For all q = (q^, q/) E W, Vp = 
(vp^jVpj,) we dehne 


w{vp,q) 



■ qm + (ni 5 ^U 25 ™)divq,„jdx 
(XpfVp^ ■ q/ + (ni?^u®^)div^q/jdr(x) 


aSx 


Inpflmi'n-PfVp^ - llVfVp^ - l^p^Vp^ 


dr(x) 


and Wp{q) = maxo^„j,gx^ lkp||p q)l- ^ sequence (T’^)ieN of gradient discretizations is 

said to be limit conforming, if for all q = (q^, q/) G W holds 


lim Wpi (q) = 0. 

l^CSD 
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Lemma 3.1 Let (r>0«eN = (X°,, and 

(^^^)zeN = nx)i^, (n^i be two sequences of gradient discreti¬ 

sations of (5) and let us assume that (Il^)zeiM is coercive, consistent and limit conforming. Let 
us furthermore assume that the sequence (C^, ^i)zgtisi, defined by 

+ \\^v)Vvij - L^(T) + X] j, 

aSx 

satisfies 

>i“Ci,.,S-=0 (6) 

and that there is a constant C G K independent of I such that 

Yl <C - YijylVjyl^\\L2(^^) (7) 

oSx «6x 

for all Vx>i G X^i, / G N. Then {V^) zeN is coercive, consistent and limit conforming. 

Proof Coercivity: is coercive, since for all / G N we have (with T> = T>\T> = T>^) 

\\LIv„,vv,^\\l‘^{q.) + \\U.VjVVf\\L‘^{T) < {Cvp + Cv)\\vv\\v < max(l, C)(Cx),© + Cv)\\vv\\v 
and since niax(l, C) ■ + Cx>i) is uniformly bounded. In the last inequality we have used 

that ||nx)||D < max(l, C)||vx)||^, which follows from (7). 

Consistency: Let / G N be hxed and V = V^V = v\ We first choose, for a given u = 

{um,Uf) G a w2? ^ ^2?) such that SD{vjy,u) = Sd{u). Using the inequality 

s^{vv,u) < sv{vv,u) + ||nx)^r;x)„ - Ilv,^vvrr,\\L2{n) + \\IlvfVVf - ^VfVVfWi'^ir) 

+ WUvfVVf - + Y, W'^Vm’^Vrr, “ (F , 

aSx 

which holds for all vv G X®, we obtain 


Sxiir) < Sx>{u) + Cx),x)IIlx)IIi’- 

Moreover 

IILdII© < Sx>{u) + \\XUm\\L‘^(nY + II^TM/||L 2 ('r)d-l + Y^ ll(7a'*^m — 

oSx 

which implies that H^xz'll®' is uniformly bonded and therefore S^i{u) —)■ 0 as / —)■ cx). 

Limit Conformity: Let again / G N be fixed and V = V^V = v\ For given q = (q^,q/) G IF 
and vt> G we calculate 

< w©(t;©,q) + \\Tivr,,Vv,^ - T\.Vr^VvJ\m^) ■ ||divq™||i2(f^) 

+ W^VfVVf — ^VfVx>f\\L^(T) ■ ||divT-qj||i2(r) + 'Y^{\\^T>fVVf — ^VfVVf\\L^(ra) 

oSx 

+ W^VfVVf - ^VfVVfWl^iFY + W^Vm'^'Drr, “ {T ■ ll7n,aqm||L2(r„) 

< Wv{vv,(\) + Cv^- W^vWv ■ (||divq,„||i2(f^) + ||div^q/||i2(r) + Y ll7n,aqm||L2(rc,))• 

aSx 
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Taking (7) into account, we derive 


>V©(qm, q/) < max(l, C) ■ sup < max(l, C) ■ q/) + Cvp\H\w)- 

Oj^v-p&Xj, \\Vv\\v 

Therefore W^i(qm, q/) tends to zero as I goes to infinity. □ 

Proposition 3.1 (Regularity at the Limit) Let (T’*)«eN be a coercive and limit conforming 

sequence of gradient discretizations and let Vx^O^eiM be a uniformly bounded sequence in 

X^i. Then, there exist {vm,Vf) G and a subsequence still denoted by such that 

in L'^iyiY, 

VvfV'Di XrVf in L‘^{TY~^, 

- ^aVm in L'^{Ta), for all a e y. 

Proof By dehnition of the norm of X)),, and by coercivity, 

and (n^^Mx)!^ — Yijy^Ujyi ), a G y, are uniformly bounded in (for I —)■ cx)). Therefore there 

exist Vm G T^(f2), Vf G T^(r), G G Lf(ytY, H G L^(r)‘^“^ and Ja G L^(rQ), a G y, and a 

subsequence still denoted by {vxn ,Vj)i )igiM such that 

m f' 



Vm 

in L^(f2), 

■ 

-- G 

in LY^Y, 


^Vf 

in LYT), 

X'DfV'pi^ — 

^ H 

in LYtY~\ 


n“ 

Vt>i .Ja in 


in T^(r«), for a G y. 


Using limit conformity we obtain (by letting I —)■ cx) 

/ (G ■ q,„ + n,„divq^)dx + / (77 ■ q/+ n/div^q/)dr(x) + V] / 7 n,aqm(-^Q - t'/)dr(x) 


oSx 


= 0 

( 8 ) 


for all (qm, qj) G x The statement of the proposition follows now from Lemma 2.1. 


Corollary 3.1 Let (T>’')i^^ be a sequence of gradient discretizations, assumed to be limit con¬ 
forming against regular test functions (qm,q/) ^ ^ ^Wf )^eiM be a uni¬ 
formly bounded sequence in X^i, such that <^nd are uniformly bounded in 

(for I ^ oo). Then holds the conclusion of Proposition 3.1. 


3.2 Application to (5) 


The non conforming discrete variational formulation of the model problem is dehned by: hnd 
{nv„,,uvf) G such that 


Tf 


/ ^mXvmUv,^ ■ V©„n2,„dx + / AfVvfUVf ■ Vv^vv^drYx) + J2 of - 1 

JT JTi S 

+ (1 - - fiv^vv^'^drix) 

(a,/3)e{(a±(i),«=F(i))} 

- [ h^IIi^^nT^^dx - f hfUjy^vVfdTf{x) = 0, 


( 9 ) 


12 



for all {vv^,vvf) e X^. 

Proposition 3.2 Let ^ G (|,1] and V be a gradient discretization, then (9) has a unique 
solution {ux>„^,ux>f) G satisfying the a priori estimate 

\\{uv,^,UVf)\\v < C[\\hm\\L^{Q.) + ||%||L2(r)j 

with C depending only on f, Cx>, A^, A/? df, df, and Aj,^. 

Proof The Lax-Milgram Theorem applies, which ensures this result. □ 


The main theoretical result for gradient schemes is stated by the following proposition: 


Proposition 3.3 (Error Estimate) Let u = {um,Uf) G , q = (qm,q/) he the solution 
of (4). Let f E {^,1], V be a gradient discretization and uv = G Xlf, he the solution 

of (9). Then, there exists Co > 0 depending only on f, Cv, A^, kf,^m, A/, df, df, Xf^^, and 
Xf^n such that one has the following error estimate: 

~ '*^m||L2(Q) + — M/||L2(r) + — n/||L2(r) 

+ ^ ~ 7a'*^m||L2(r„) + \\X ||+ \\XrUf — ||i2(r)d-l 

a€x 

< Co(Sv(um,Uf) + Wviqm, q/))- 


Proof From the dehnition of Wd, and using the dehnitions (4) of the solution u, q and (9) of 
the discrete solution ut>, it holds for all (vx>^,vx>f) G 


> 


E 

iei 


E 

i&I ' 


■ qm + {Tlvrr.vvjh^d-K + J^(XvfVVf ■ Q/ + {IlT)^vvf)dfhf^dT{x) 

(^^laUm + (1 - OlgUm “ «/) j dr(x) 


f _Ti_ 

Ir, 2^-1 


(Q,/3)e{(a±(i),aT(i))} 

- Vn,„))dx + f - V^M/))dr/(x) 


Tf 


2^-1 


Y - n: 




(«,/3)e{(«±(i),«T(i))} 


X U-faUm + (1 - - Uf - - (1 - J dr(x) 


0 ) 


Let ns choose tc® = {wx>^,wx>f) G s.t. s{wx>,u) = Sx>{u) and set {vx>^,vx>f) = ux> — u)x> in 
(10). Then holds 


llVUm — + ||VrM/ — ||i2(r)d-l 

+ Y - ^VfUVf - IcUm + M/||L 2 (r^) < C ■ {Sv{u,^, Uf) + >V©(qm, q/)), 

aSx 

with a constant C > 0 depending only on Am; kfAm, A/, df, df, A/,„, and A/^„. Taking 
coercivity into account leads to the statement of the proposition. □ 
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4 Two Examples of Gradient Schemes 

Following [7], we consider generalised polyhedral meshes of Let M be the set of cells that 
are disjoint open subsets of such that iF = For all K G M, denotes the 

so-called “center” of the cell K under the assumption that K is star-shaped with respect to 
^K- Let denote the set of faces of the mesh. The faces are not assumed to be planar for the 
VAG discretization, hence the term “generalised polyhedral cells”, but they need to be planar 
for the HFV discretization. We denote by V the set of vertices of the mesh. Let Vx, Vo- 
respectively denote the set of the vertices of iF G Ad, faces of iF, and vertices of a G For 
any face a G we have Vo- C Vk- Let (resp. J^s) denote the set of the cells (resp. faces) 
sharing the vertex s G V. The set of edges of the mesh is denoted by £ and £„ denotes the set 
of edges of the face a G Let Te denote the set of faces sharing the edge e G T, and Ado- 
denote the set of cells sharing the face a G A”. We denote by Text the subset of faces a G A 
such that Ado has only one element, and we set Eext = Uo-eJ'e^rt ^ The 

mesh is assumed to be conforming in the sense that for all a G A \ Text) the set Ado- contains 
exactly two cells. It is assumed that for each face a & T, there exists a so-called “center” of 
the face Xo- such that 

Xa = ^ Xs, with ^ = 1, 

•SGVct SGVct 

where /do-,s > 0 for all s G Vo-. The face a is assumed to match with the union of the triangles 
To-,e dehned by the face center Xo- and each of its edge e G To-. 


The mesh is assumed to be conforming w.r.t. the fracture network F in the sense that there 
exist subsets Ar^, i G / of A such that 

r.= U if. (11) 

(TeJ'r- 

We will denote by Ar the set of fracture faces Ar^. Similarly, we will denote by Tr the set 
of fracture edges IJo-eJ'r and by Vr the set of fracture vertices IJo-eJ'r 


We also dehne a submesh T of tetrahedra, where each tetrahedron DK,cr,e is the convex hull 
of the cell center x^^- of K, the face center Xo- of a G Tk and the edge e E £a- Similarly we 
dehne a triangulation A of F, such that we have: 

A = IJ DK,a,e and A = (J A,e. 

We introduce for A G T the diameter ho oi D and set hj- = maxoer ho- The regularity of our 
polyhedral mesh will be measured by the shape regularity of the tetrahedral submesh dehned 
by dj- = max/)g 7 - — where p_D is the insphere diameter of A G T. 


The set of matrix x fracture degrees of freedom is denoted by x dofjy^. The real 

vector spaces and of discrete unknowns in the matrix and in the fracture network 
respectively are then dehned by 


where 


Xv^ = span{e^ | z/ G dofjy^} 
Xvf = spanje,, | u G dofj)^}, 




for z/ G dofjy^ 
ioTuedofjy^. 
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For uj)^ G Xj)^ and v G we denote by Ujj the z/th component of ux>^ and likewise for 

G Xdj and u G dofjy^. We also introduce the product of these vector spaces 

Xj) = Xj)^ X Xtdj, 

for which we have dimXx) = + 4^dofjy^. 

To account for our homogeneous boundary conditions on (9^2 and Sq we introduce the subsets 
dofoirm C dofjy^, and dof C dofj^^, and we set dof^ir = dof x dof^irf, and 

= {m G Xd \ u^ = 0 for all u G dof 

4.1 Vertex Approximate Gradient Discretization 

In this subsection, the VAG discretization introduced in [7] for diffusive problems on hetero¬ 
geneous anisotropic media is extended to the hybrid dimensional model. We consider the Pi 
finite element construction as well as a hnite volume version using lumping both for the source 
terms and the matrix fracture fluxes. 

We hrst establish an equivalence relation on each At*, s G V, by 

K =Ms there exists n G N and a sequence in Xs\Xr, 

such that K G Mai-,L G A1<t„ and Ado-i+i H Mat 7^ 0 
for i = 1,..., n — 1. 

Let us then denote by Xig the set of all classes of equivalence of Xis and by Ks the element 
of Ais containing X G Ad. Obviously Ais might have more than one element only if s G Vp. 
Then we dehne (cf. hgure 2) 

dofjy^ =AiU I a G Ar, a: G Ai^j U |x, | s G G Af,|, 

dofvf = At U Vr, 

d(^f Dirm ' 'Ia's I S G "Vexti Kg G , 

dof Vr n Vext- 

We thus have 

= ^Uk I X G Ad| U I <J G Ar,K G Ado-j 

U s G V,X, G Ad,|, (12) 

Xvf = jwa I (y e Jt| U I s G Vrj. 

Now we can introduce the piecewise affine interpolators (or reconstruction operators) 

Hr: ^ (0\r) and Ha ; Xj,^ H\T), 

which act linearly on Xjy^ and Xx)^, such that IlrUVm is affine on each DK,a,e ^ and satisfies 
on each cell X G Ad 

= Ur, 

nrM©^(xs) = Mjr, Vs G Vir, 

nr'*^»m(xo-) = ur^ V(t g Ar n Ar, 

^rUVrr^M = E l3a,sUR^ Va G Ar\At, 

S&Va 
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Figure 2: Cell K touching a 
fracture face a. Illustration 
of the simplices on which: 
Red: is constant. 

Grey: is constant. 



while is affine on each T^-^e ^ and satishes for all v G do/x)^ 

n^MX)j.(Xjy) Mj/, 

where G is the grid point associated with the degree of freedom v G dofx>^ U dofj)^. The 
discrete gradients on Xx>^ and Xti^ are subsequently dehned by 

= VIIt- and = VtIIa. (13) 

We dehne the VAG-FE scheme’s reconstruction operators by 

• = nr, 

• = Ra, (14) 

• = 7„nr for all a G y. 

For the family of VAG-CV schemes, reconstruction operators are piecewise constant. We 
introduce, for any given K E AA, a. partition 

K = ukg(^ U ( U 

s£VK\Vext crer/fnrr 

Similarly, we dehne for any given a G Jr a partition 

a = u„u(^ IJ j. 

SGVcrVVext 

With each s G V \ Vext and Kg G AAg we associate an open set uJks^ satisfying 

^Ks = U 

KeKs 

Similarly, for all s G Vr \ Vext we dehne Ug by 
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We obtain the partitions 

= (^ IJ , r = (^ IJ u^y 

Uir^ !^edo/i,^\do/oir^ 

We also introduce for each T = T„^s,s' € A a partition T = which we need for 

the dehnition of the VAG-CV matrix-fracture interaction operators. We assume that holds 
1^11 = 1 ^ 2 ! = l^sl = ||T| in order to preserve the hrst order convergence of the scheme. 

Finally, we need a mapping between the degrees of freedom of the matrix domain, which 
are situated on one side of the fracture network, and the set of indices y. For K^j G dofjy^ 
we have the one-element set x(A'o-) = {o; G x | = Hq on a} and therefore the notation 

a{K^) = a G x(A^). 

The VAG-CV scheme’s reconstruction operators are 

• 

• Ilx>fUT>f = Uutuj,,, 

' ' (15) 

• Uj^fUDf = ^ {UctIti + UstT2 + Us'Its), 

G.s.vsA 

• = X] X] {uK,lT^+Uj^^lT2+U^^,^T3)daiKa)a^ro,- 

Remark 4.1 The VAG-CV scheme leads us to recover fluxes for the matrix-fracture interac¬ 
tions involving degrees of freedom located at the same physical point (see subsection 4-3)- 

Proposition 4.1 Let us consider a sequence of meshes (Al^)zeiN and let us assume that the 
sequence (T^);gtisi of tetrahedral submeshes is shape regular, i.e. Oq-i is uniformly bounded. We 
also assume that hm/_,.oo hq-i = 0. Then, the corresponding sequence of gradient discretizations 
defined by (12), (13), (14), is coercive, consistent and limit conforming. 

Proof The VAG-FE scheme’s reconstruction operators are conforming, i.e. Vt> C V°. There¬ 
fore we deduce coercivity from Proposition 2.1. Furthermore we have by partial integration 
>V©(qm,q/) = 0 for all (q^, q/) e W. Hence (V‘) ;giM is limit conforming. 

To prove consistency, we need the following prerequisites. We define the linear mapping 
such that for all fjm G Off and any cell K ^ A4 one has 

{PVm''P'm) K — f^raip^K)-, 

(^©„^m)x^ = ^m(x,) Vs G Vk, 

(,PT>m^rn}Kcr V(J G XK G VV’ 

Likewise, we dehne the linear mapping Pvf- —)■ such that for all G holds 

{Px>jf)f)y = ijjf(xi,y) for all u G dofvf It follows from the classical Finite Element approximation 
theory and from the fact that the interpolation fla,s{Pvrr,'4’m)K^ at the point x^-, a G pK\Pr 

s&Va 

is exact on cellwise affine functions, that for all x holds 

W'P.TPvrr.i’m - '«/’m||//i(n\r) + W'^APvffl’f - fl’fWmir) < Cfipm-, fl’f, 0r)hr- (16) 
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The trace inequality implies that for all v G Hl^{yL\r) holds 

\\lav\\L^ir^) < C'(ff\r)i|n||^i(f2\r) for « ^ y. 

We can then calculate for {um,Uf) G x 

<Sv{um-,Uf) < y/2\\IlrPv^Um - ■Wm||iyl(n\r) + WlaO^TPVmUm “ Mm)||L2(r„) 

«ex 

+ X] Vs\\IlAPvfUf - Uf\\m{ri) 
i£l 

< C{n\T, #x, #/, {Um, Uf), Or) hr- 

Since x is dense in the sequence of VAG-FE discretisations (T’^)igtisi is consistent 
if h'fi —)■ 0 and 9-^-1 is bounded for / —)■ cx). □ 

Proposition 4.2 Let us consider a sequence of meshes (Wl^)zgiN and let us assume that the 
sequence (T^)zeiisi of tetrahedral submeshes is shape regular, i.e. 6q-i is uniformly bounded. We 
also assume that hm;_^oo hq-i = 0. Then, any corresponding sequence of gradient discretizations 
defined by (12), (13), (15), is coercive, consistent and limit conforming. 

Proof We combine Lemma 3.1 and Proposition 4.1. Thus, we have to show that the assump¬ 
tions of Lemma 3.1 are satished, where (V );giM corresponds to the sequence of VAG-CV gradient 
discretisations and (T’^)zgtisi to the corresponding sequence of VAG-FE gradient discretisations. 

For the following, we dehne = U.«. and V“ = Va- To ease the notation in the 

proof, we will use, for a G y, the uniquely identihed mapping /i“: V" U C dofj^^ — )■ dofjy^, 
dehned by pP{a) = (such that x{Po-) = {«}) and /i“(s) = Kg (for a cell K such that 

K G AAa with (j G n Ps and x(i^o-) = {«}). Let now a G y be hxed. Since the mesh is 
conforming with respect to the fracture network, there is for every a G e = ss' E a, 
z/((T, e) G {(T, s,s'}, such that 

sup - IlVfVVf){^)\ = - flvWVf){^u{a,e))\ = “ Vuia,e)\- 

XST^.e 

Then we have 


< 'y ^ y \Ta,e\\ViJ.°‘(u{<T,e)) n,y(o-,e) | 

We have to check (6) now. It can be verihed that [4], Lemma 3.4 applies to our case, both, in the 
matrix domain, where face unknowns might occur, as well as in the fracture network, a domain 
of codimension 1. This means that we can state that there exist constants CmiO-r), Cf{6-j-) > 0, 
such that 


\\P-v,^uv^ - Vv^vv,^\\L-^(n) <Cm- hr ■ \\'^v,n'>^Vm\\L2(nA und (17) 

W^VfVVf — ^VfVVf\\L2{r) = W'P-T'fVrf — 'AvfVVf\\L2{T) < Cj ■ La ■ IIVD^WD^||p^2(r)d-i (18) 

For the following calculation we take into account [4], Lemmata 3.2 and 3.4. We also use that 
the mesh is conforming with respect to the fracture network and that for a E P and K E 
(or equivalently for K E M, a E Pk) holds: hx is asymptotically equivalent to h^ and |iF| is 
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asymptotically equivalent to ha\a\^ where hx '■= maxT-^ocic^D and := maxAgTco-^r- Let 
a G X, a G and K G such that x(iLo-) = {a}. Then we have 

= II ^ - n®^e,,)||^2(^) 

!^e{cr}U(VCT) 

— L* ■ 1*^1 ^ ^ ~ ‘V^°‘{a)) 

s&Va 

< L* ■ ■ f - vkY + < c* ■ ■ ||V©„n©„||^2^^^d. 

^ sGVx aG^^n^r 

Therefore 

- ^V^'VvJ\l2^^„) <C -hA- \\^VmVvJ\l2^^f- 

(19) 

Altogether we obtain 

l|LlD„t^©™ - nD„nD„||L2(Q) + \\YlvjVvj - nD^nD^||L2(r) + WUtdjWj - nD^nD^||L2(r) 

+ Yj ~ Llp^nx)™||L2(r^) < C ■ {hr + hA + h^) ■ \\{vT>m,'^Vf)\\v, 

aSx 

with a constant C depending only on and Or- This proves that (6) is satished. □ 

1 

Corollary 4.1 The precedent proof shows that Sx>{um,Uf) = 0{hf) for {um,Uf) G x C“ 

1 

and that Wri^im, Q/) = 0{hf-) for (q^, q/) G x . However, we can prove a higher order 
of convergence, i.e. W©(qm,q/) = 0{hr) for (qm,q/) G x and Sr{um,Uf) = 0{hr) 
for {um,Uf) G X C^. 

Proof Consistency: Classically, for all {iprn,Tf) ^ x C“, we have the estimate 

W^Vr^PVr^Tm - V^m||L2(0) + l|Ll®^pD„V^m - 7o<dm||L2(r^) 

+ \\P-VfPvfTf ~ Tf\\L^{r) + WP-VfPvfTf ~ 9^/IU2(r) < cst{ipm, Tf) ' hr, 
while (16) grants that holds 

\\'^VruPvrr,Tm ” V111,2+ \\VvfPvfTf ~ V99||L2(r) < cst{(pm,Tf,hT)hr- 

Taking into account that Cff x is dense in V, we see that the treated discretisation is 
consistent with Sx>{Tm,Tf) = ^{hr) for {iprn,Tf) ^ ^ 

Limit Conformity: For all T G A and for all G Ax)„, we have that 

J^{Ilf,^ur,^ - Ul^uvjdrix) = 0. 

Introducing the linear operator P : L^{Ta) —t L‘^{Ta) such that Pi’p) = /^y:>dr(x) on T for 

all T G A, we hrst calculate for any q^ G C^^ 

hn.aqm “ ^(hw^qm) ||i2(r^) = Y Y IlFw^qm - ^(7n,aqm)||i2(r) < C{qm,0r) ' h^- 

(t^Tol AbTGct 
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We proceed: 


1/ 7n,aqm(n^„M©^ - n^^M7,„)dr(x)| 

JTc, 

= I / (7n,aqm “ P(7n,aqm)) - n^„M25^)dr(x) | 

^ ||7n,aqm ~ P{'ln,a^.m)\\L^{Ta)\\^'DTri^'Drri ~ ^'Dm^T^m\\L^{^a) 

3 

< (^(qm, 

for all qrn ^ where we have used (19) in the last inequality. We can now conclude by 

calculating for all for q = (q^, q/) G x C(y^ 

w^{uv, q) = (w® - wv){uv, q) 

= / divqm(n©„ - n©„)u®„dx + / div^q/(n®j. - n©^)M©^dr(x) 

JO Jr 

+ ^ / 7n,aqm(^(nD^ - nx)^)-Ux)^ - (Uj)^ - nx)^)-Ux)^ - - n^^)-Ux)^jdr(x) 

< l|nD„Mx)„ - nx)„-ux)^||i2(j^) ■ ||divqm||2,2(Q) 

+ W^VfUVf — ^VfUVfWi^ir) ■ ||div^q/||L2(r) + f(UllpfUpf — ||i,2(r^) 

aSx 

+ \\P-Vf'U'D^. — ||L2(rc«)) ■ Il7n ,a^m lU^r.) 

+ [ 7n,aqm(n^^M©„ - U'^^uvJdT{x)] < CiOr, q) ' ' \\uv\\v, 


where we have taken into account the conformity of V in the hrst equation and (17), (18) in 
the last inequality. □ 


Remark 4.2 The proofs of Propositions f.l and 4-2 show that for solutions {um,Uf) G and 
(qm,q/) eW of (4) such that Um e C^iK), Uf G C‘^{a), q^ e {C^{K)Y, q/ G {C^{a)Y~^ for 
all K & M. and all a E Tf, the VAG schemes are consistent and limit conforming of order 1, 
and therefore convergent of order 1. 


4.2 Hybrid Finite Volume Discretization 

In this subsection, the HFV scheme introduced in [8] is extended to the hybrid dimensional 
Darcy flow model. We assume here that the faces are planar and that x^. is the barycenter of 
a for all a E P. 

The set of indices dof^y^ x dofjy^ for the unknowns is defined by (cf. figure 3) 

(io/i,„ = XI u(|jA 7 ,) 

cfGlT 

dofvi = Pr A Tr, 

Dirm ~ -Pe^t i 
dof Tp n Text) 
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Figure 3: Cell K touching a 
fracture face a. Illustration 
of the polyhedron and poly- 
gone on which: 

Red: is constant. 

Grey: is constant. 



where ioi a E and K E 


f M. ifaEjr\J^r 

" \ {K} ifaEJ^r- 

and = {K^ \ K E Ada}- We thus have 

I ^ ^ Wl| U I cr E Ar,Ka ^ Wlcrj, 
Xx>j: = 1^(7 I cr E -^rl U | e G £^r|- 


( 20 ) 


The discrete gradients in the matrix (respectively in the fracture domain) are dehned in 
each cell (respectively in each face) by the 3D (respectively 2D) discrete gradients 


(resp. Vdj.) as proposed in [8], pp. 8-9. (21) 

The function reconstruction operators are piecewise constant on a partition of the cells and 
of the fracture faces. 

These partitions are respectively denoted, for all iF G Al, by 

K = UJk U ^ 


and, for all cr G Ap, by 

a = U IJ 07,7,e)- 

\^erct 

With each cr G A \ Aext and iF^- G Alo- we associate an open set s.t. 

= U 

KeK^ 

Similarly, for all e G \ ^ext we dehne cOe by 
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We obtain the partitions j T = ) • 

We also need a mapping between the degrees of freedom of the matrix domain, which 
are situated on one side of the fracture network, and the set of indices x- Foi’ cr E J^r and 
Ka G Aia holds by dehnition = {K} for a K E Ma and hence = nx,<T is well dehned. 
We obtain the one-element set xi^a) = {a E x \ = rio and therefore the notation 

a{K^) = a E xi^a)- 

We dehne the HFV scheme’s reconstruction operators by 


uedofj,^\dof 

• Hx’f'U'Dj = 

u&dofj,^\dof 

Dir j: 

• K(K^)aUK^K for all a e X. 


( 22 ) 


Proposition 4.3 Let us consider a sequence of meshes (Wl^);giH and let us assume that the 
sequence (T^);gtisi of tetrahedral submeshes is shape regular, i.e. Oq-i is uniformly bounded. We 
also assume that hm;_^oo hq-i = 0. Then, any corresponding sequence of gradient discretizations 
(P^)igiM, defined by (20), (21) and definition (22), is coercive, consistent and limit conforming. 

Proof Let us denote in the following by fix and Iljr = Iljr the HFV matrix and fracture 
reconstruction operators for the special case that = 0 = cUe for all E Uo-eJ'*^^ 
e G £^r- We start our numerical analysis for HFV by proving the proposition for these special 
choices and then use Lemma 3.1 for generalizing the results. 

Coercivity: We hrst prove that limit conformity against regular test functions, as proved 

below, implies coercivity. 

Assume that the sequence of discretizations (F’*);giH is not coercive. Then we can hud a 
sequence {{ux>i ,Ujyi ))«gN with ,Ujyi ) G Wll,, such that 

m m U 

\\'nvi,,UjyiJ\L^^n) + \\Uj,iu^\\L^(^r) = 1 and \\{ut) 1 ^,u^)\\vi < y. (23) 

Then follows from a compactness result of [21] that there exists a m = {um,Uf) E L^(H) x L^(F), 
s.t. up to a subsequence 

^ (wm, Uf) in L‘^{n) x ^^(F) ( for / ^ oo) 

and therefore ||Mm||L2(o) + ||'*^/||L 2 (r) = 1. On the other hand follows from the discretizations’ 
limit conformity against regular test functions (see below) by Proposition 3.1 and Corollary 3.1 
that {um,Uf) E V° and that up to a subsequence 

f Vn™ in L‘^{nY, 

J mL\rY-\ 

[ ^ w - in ^^(F^), for a G x- 
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Since by construction holds \\{uj)i^,U'pi^)\\j)i —)■ 0, we obtain \\{um,Uf)\\vo = 0. But || ■ ||yo is a 
norm on which contradicts the fact that ||Mm||L 2 (o) + ||'W/||L 2 (r) = 1. 

Consistency: For {(pm,^f) G x let us dehne the projection Pv^ipm £ 
for all cell K & Ai one has 

{Pv^^m)K = <^m(xx), 

{P'Dm^rnjK^ — VcT E. PK■> 

and the projection Pvf^^f ^ such that (Px)/V?/)y = ^f(Pu) for all z/ G dofx>f Let us set 
Vv = {Pvm^m,Pvfy^f)- Then holds 

\\vk - ‘fm\\L^(K) < ■ hr ■ \K\^ for K E M, 

where := max^ ||V99m||. Summing over K E M. yields 

\\P-mVv^ - 9^m||L2(o) < ■ hr ■ |hl|2. 

We also have 

\\vk^ - 7 a<^m|U 2 (r,) < c“^ ■ hr ■ |ct| 5 for a G P", E Af“ 
where c“^ := maxp^ || VT-7Q,93m||, from which we obtain 

- 7a<^m||L2(r,) < ' hr ■ |r„|i 

Analogously we can derive 

\\P-TVrj - V^/||L2(r) < c^j ■ hr ■ |r|2, 

where := maxp ||Vr9?/||. Furthermore, it follows from Lemma 4.3 of [8] that there exists 
C > 0 depending only on 9r and (p such that 

II - V(p\\l‘ 2 {p,) + II - V99||L2(r) < Chr 

Taking into account that x is dense in V^, we see that the treated discretisation is 

consistent. 

Limit Conformity: Let (pm ^ ^Wm for all iF G Ad, a E Pk let Pk '■= <Fmdx and 

PK,a '■= Xj 7 n^<^<Fmdr(x). lu exactly the same manner as [19], (29)-(31) are proved, we can 
show that holds 


^vJJ^'Drr, <Chr\\VvrnUv^\\L 2 (^r)d and (24) 

«6x 

= 5^ ^\(^\{uK-UK^){pK,a-PK)-nK,a, (25) 

K&M ctST'k 
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where 


KeM aeJ^K 

= X] X] RK,a{uVrr,)^K,a ' / <^mdx, 
K&MaeTx JDk„ 

= E E \<^\uK'^K,cr ■ ^K,c 

KeM a&J^K 






K ,a 

and 


= f ■ (fim + OlMUvJdiv{(p^)^dx, 


with the dehnition of the gradient stabilization term RK,a{uvm) in [ 8 ], pp. 8-9. Therefore, 
applying Canchy-Schwarz ineqnality to (25), using the regularity of cpm, and the estimate (24), 
we deduce that there exists C depending only on cpm, such that 


■ (Pm + (JlMUvJdw{(pr^ 


a&X ^ 


7n,a‘^m(n^„M©™)dr(x) < Chr\\Vv^uv„ 


Taking into account the result [19] (33), i.e. for all (p G exists a constant C > 0 depending 
only on 6 ^ 7 -, such that 


we obtain all together 


'^VfUVf ■ Pf + (nj-M®^)div(c^/)jdr(x) 

E Ch/\\\Wx>jyUx> J:\\^ 


wv{uv, q)<C-hr- \\ut)\\v for all q G x 

This result is shown above to imply coercivity, which is needed to conclude now. 

Finally, using that x is dense in IF and the coercivity of the scheme, we derive 
limit conformity on the whole space of test functions. 

Generalization to arbitrary HFV discretizations: We want to apply Lemma 3.1. From [ 8 ] 
Lemma 4.1 and [21], it follows that there are positive constants Cm and Cj only depending on 
9-j- and d, such that for all ux> G Xjy holds 

W^M^Vr^ - nD„-Ux)„||i 2 (Q) = E E — R'm ■ hr ■ ll^»m'*^n’mlli 2 (Q)d 

K&M (t^Tk 

Ilhlj'WD^ — ||L 2 (r) = \oj^^e\{Ua- — U^f ^ Cf ■ h\ ■ || || ^ 2 (r)‘*-l' 

The remaining conditions of Lemma 3.1 are trivially satished, from what follows the statement 
of the proposition. □ 

Remark 4.3 The precedent proof shows that for solutions {um,Uf) G and (qm,q/) ^ hF 
of (4) such that Um G C‘^{K), Uf G C‘^{a), q^ G {C^{K)Y, qj G {C^{a)Y~^ for all K ^ Ai 
and all a eT f, the HFV schemes are consistent and limit conforming of order 1, and therefore 
convergent of order 1. 


L^(nf- 
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4.3 Finite Volume Formulation for VAG and HFV Schemes 

For iF G At let 


^ f — ) s G Vx} U {iFo-, (j G AA n Ar} for VAG, 
- S for HFV. 


Analogously, in the fracture domain, for a G Ap let 

do fa = 

Then, for any z/ G dofx the discrete matra-matra-fluxes are dehned as 


Va for VAG, 
Sa for HFV. 


HKuiH'Dm) ^ ^ i I AyyjVx)„^C;/Vx)^Ci/'dx j ("Utp tiy'). 

u'^dof^^dK ^ 


such that = Y.k&m '^u&dofK Fku{uvJ{vk - v^). For all z/ G dofa 

the discrete fracture-fracture-Huxes are dehned as 


FauiuVf) — y AfV'Pj.c^V'p^Cjy'dTflj 

u'edou 


X [Ua - u,/, 


such that J^Af'VvfUVf'VvfVVfdTjix) = Yl,a&Fr'^y&dof^ Fai,{uv^){ya-Vi,). To take interactions 
of the matrix and the fracture domain into account we introduce the set of matrix-fracture {mf) 
connectivities 

c = {(z/^, Vf) I Vm e dofjy^, Uf G dofj)^ s.t. x^^ = x^J 

with dof^^ = {z/ G dofjy^ | x,^ G F}. The m/-huxes are built such that 

('i'XZ™ Axz/)) = ^ F,,^yj{U'D^,U'Dj){Vy^-Vtyj) 

= T.f ^ E + (1 - - ni,,«i,,)dr(x), 

i&I d'^i ^ {a,l3)£ 

{fo±(i),aT(i))} 


for all G Xx>. For all a G Ap and K G Alo-, let us denote by a{K, a) the unique a G y 

such that a & Fa and Hq, = nx,o-. Let us also set for all a G Ap, (y x x)a = {{c({K, a), a(L, a)), 
(a(L, a), a(K, a))} with = {K, L}. Then, holds 

E E / ^hny “®,. + (1 - - nv,uv,) (ny «!,„ - ni,,«i,,)dr(x)^ 

aeJ-r (a,/3)e(xxx). ^ 

For all a G Ap, K G Ai^ and x G a, let us notice that, for the VAG scheme, one has 
Cip^(x) = lijyjtaix), and = nx)^es(x) for all s G Vo-, and for the HFV 

scheme, one has n^^^’'^^e;i^^(x) = nx)^eo-(x) = It result after some computations that the 
VAG matrix fracture huxes are dehned by 

FK^a{uv,^,UVf) = ( / x^^(nfoeo)(nfoe,)dr(x)) + (1 - “ fo) 

seVa ^d^ ^ 

+ (/ ^^^(nfoe^) 2 dr(x)) + (1 - Oul, - , 
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for all a G Jr, -Ma = {K, L} , and by 


FqJuv^,uv,)= y1 J2 { 

f^edJcjgQ^ J'Q)nJ'.nJ'r KeM^nQ^,LeMA{K} 

(j + (1 - O^L. - Us^ 

+ Y1 if + (1 - 0^1^, - Ms') 

s'ev,, I ss'££^ ^ 

+ (/ ^^^(nn^ ea)(nc^e,)dr(x)) + (1 - }, 

for all s G Vr, Qs G A'ls- Similarly the HFV matrix fracture fluxes are defined by 
F-K^AuVrn^uv^) = T’/Wdr(x)) + (1 - -m^), 

for all a G Jr, = {K, L}. 

We observe that for the VAG-CV scheme (since T/-(nx)^es/)(nx)^.Cs)dr(x) = 0 for s 7 ^ s' 
and T/(nx)^eo-)(nx)^es)dr(x) = 0) as well as for the HFV scheme, the fluxes only 

involves the d.o.f. located at the point Xj,^ = Xy^.. 

The discrete source terms are dehned by 



Figure 4: mm-fluxes (red), m/-fluxes (dark red) and j^-fluxes (black) for VAG (left) and HFV 
(right) on a 3D cell touching a fracture 


The following Finite Volume formulation of (5) is equivalent to the discrete variational 
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formulation (9): find {ux>^,Ux>j) G such that 

^ foralli^eAf; E Fku{uvJ = Hk 

u&dofK 

for all u G Jt : E F„u{uvj) - E F^^„{uv^,uvj) = 

U&dofa t'm&dof-p^ 

S.t. {Um,0-)£C 

for all z/„ G dofv^ \ (A^ U dof : 

- E FKu^{UVru)+ Y. F^^^^.{UVm,UVf) = 

K&Mv^ UfGdofTij: 

s.t. {Vm.,Vf)&C 

for all Vf G dofvf \ {Fr U dof : 

~ E Fa;/^{U'D^) — E Fy_^y^{U'P^,U'P^) = Hyj. 

o'GJ'r , 1 / j &dofT>^ 

s.t. (Vm,Vf)&C 

Here, M.v^ stands for the set of indices {K G A1 | z/m G dofx} and Fr,vf stands for the set 
{a e Fr \ i^f e dof^}. 

It is important to note that, using the equation in each cell, the cell unknowns uk-, K ^ Fi, 
can be eliminated without £ll-in. 


5 Numerical Results 

The objective of this numerical section is to compare the VAG-FE, VAG-CV, and the HFV 
schemes in terms of accuracy and GPU efficiency for both Gartesian and tetrahedral meshes on 
heterogeneous isotropic and anisotropic media. For that purpose a family of analytical solutions 
is built for the fixed value of the parameter = 1. We refer to [12], [3], [2] for a comparison 
of the solutions obtained with different values of the parameter ^ G [|,1] with the solution 
obtained with a 3D representation of the fractures. 

Table 1 exhibits for the Gartesian and tetrahedral meshes, as well as for both the VAG 
and HFV schemes, the number of degrees of freedom (Nb dof), the number of d.o.f. after 
elimination of the cell and Dirichlet unknowns (nb dof eh), and the number of nonzero element 
in the linear system after elimination without any fill-in of the cell unknowns (Nb Jac). 

In all test cases, the linear system obtained after elimination of the cell unknowns is solved 
using the GMRes iterative solver with the stopping criteria 10“^°. The GMRes solver is pre¬ 
conditioned by ILUT [25], [26] using the thresholding parameter 10“^ chosen small enough in 
such a way that all the linear systems can be solved for both schemes and for all meshes. In 
tables 2 and 3, we report the number of GMRes iterations Iter and the GPU time taking into 
account the elimination of the cell unknowns, the ILUT factorization, the GMRes iterations, 
and the computation of the cell values. 

We ran the program on a 2,6 GHz Intel Gore i5 processor with 8 GB 1600 MHz DDR3 
memory. 

5.1 A class of analytical solutions 

We consider a 3-dimensional open, bounded, simply connected domain D = (—0.5, 0.5)^ with 
four intersecting fractures P 12 = {{x,y,z) G D | a; = 0, y > 0}, P 23 = {{x,y,z) G D | y = 0,a; > 
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0}, r34 = {{x^y^z) G f 2 I a; = 0, 1/ < 0} and = {{x,y,z) G f 2 | y = 0, a; < 0}. We also 
introduce the piecewise disjoint, connex subspaces of = {{x,y,z) G | y > 0, a; < 0}, 

^2 = {{.x,y,z) G I y > 0, a; > 0}, fls = {{x,y,z) G | y < 0, a; > 0} and = {{x,y,z) G 
f2 I 1 / < 0, a; < 0}. 


Derivation: For {um,Uf) G V, we denote Um{x,y,z) = Ui{x,y, z) on Hi, i = 1,A and 
Uf{x, y, z) = Uijiy, z) on F^, ij G J, where we have introduced J = {12, 23, 34,14}. We assume 
that a solution of the discontinuous pressure model writes in the fracture network Uij{y,z) = 
af{z) + /3ij{z)'jij{y), ij G J and in the matrix domain 

{ui{x,y,z) = ai{z)ui 2 {y,z)ui 4 {x,z) 

U2{x,y,z) = a2{z)ui2{y,z)u2ii{x,z) 

Uii{x,y,z) = aii{z)uiii{y,z)u 2 ii{x,z) 

UA{x,y,z) = ai{z)u^i{y,z)uiA{x,z). 


On G J we assume 7jj(0) = 0, such that the continuity of Uf is well established at the 

fracture-fracture intersection, as well as 7ij(0) = 1, to ease the following calculations. For 
i = 1,..., 4 let iFj = and for ij G J let Tij = From the conditions 7n,aqm = 

Tf{-yaUm — Uf) on Fq,, a G x, we then get, after some effort in computation. 


«i(^) = [af{z) 
asiz) = (af{z) 
Mz) 


K 






I^l2{z, 


14 


-1 


I3ii{z) = 


KlyK2,yK4^xTl2 

KlxK^xKiyT23 
KlyTl2 

KixTi4 


KiyK^yK4xTi2 
KlxK4yT23,Ts4 


I3i2{z] 


-1 


I3i2{z), 


-/3i2{z), 


a2(z) = (^af(z) 
a4{z) = (^Q:f{z) 

(^^iz) = 

K\yK2xK'4yK4x 
^lx^2y^3x^4y 


KlyK4xTi2 


= 1 . 


KlyK2xK‘4yK4x 

KlxK^xK4yT23 
KlyK4xTl2 

KlxTi4T^4 


f3l2{z] 


I^l2{z, 


I^l2{z) 


-1 


(26) 

Obviously, we have taken af and /di 2 as degrees of freedom, here. However, these functions 
must be chosen in such a way that ^ 0 for i = 1,..., 4. 


Remark 5.1 We would like to explicitly calculate the jump at the matrix-fracture interfaces 
for this class of solutions. At Tij we have 


Uiif), y, z) - uj(0, y, z) = {ai{z) - aj{z)) ■ af(z) ■ Uij{y, z), for ij G (12, 34} 

Ui{x, 0, z) — Uj{x, 0, z) = {ai{z) — aj(z)) ■ af{z) ■ Uij{x, z), for ij G (23,14}. 

From (26), we observe, that the pressure becomes continuous at the matrix-fracture interfaces, 
as the Tij tend to oo uniformly. 


Remark 5.2 In order to obtain solutions with discontinuities at the matrix-fracture interfaces, 
we had to omit the constraint of flux conservation at fracture-fracture intersections. 


5.2 Test Case 

We dehne a solution by setting af{z) = fli 2 {z) = —1, 712 ( 1 /) = cos{27ry) + y — 1, 

'l 23 {x) = X, 734 ( 1 /) = —e“^02/) 714 ( 2 ^) = parameters we used for the different 

test cases are 



Isotropic Heterogeneous Permeability: 


Klx = Kiy = Kiz = 1 , K2x = K2y = K2z = 100 , 

Kzx = Ksy = Ksz = 3, K4x = Kiy = Kiz = 40, 

Ti 2 = 1, P23 = 0.2, T34 = 100, Ti4 = 10, 

Ki 2 = 1, K 23 = 2, = 3, i^i4 = 10. 


• Anisotropic Heterogeneous Permeability: 

Kix = Kiz = 1, Kiy = 50, K 2 X = K 2 Z = 2, A'2y = 100, 
K'iy = K^z = 3, K^x = 30, Kiz = 4, Kix = Kiy = 40, 
Ti2 = T23 = T34 = Ti4 = 1, 

7^12 = K 23 = K^i = Kii = 1 . 


In the following figures we plot the normalized L?' norms of the errors, which are calculated as 
follows: 


• normalized error of the solution: erfsoi 


\\Um Ih2(n)+||«/lh2(p) 


• normalized error of the gradient: err grad 


m 

L2(n)'i + ll^^'“/llL2(p)d-1 


In the following tables is additionally found the normalized error of the jump: 


Wg HfV 


1 Hexahedral Meshes 

Key 

Nb Cells 

Nb dof 

Nb dof el. 

Nb Jac 

Nb dof 

Nb dof el. 

Nb Jac 

1 

512 

1949 

1437 

31253 

2776 

2264 

20696 

2 

4096 

11701 

7605 

178845 

19248 

15152 

150320 

3 

32768 

79205 

46437 

1154861 

142432 

109664 

1141856 

4 

262144 

578245 

316101 

8152653 

1093824 

831680 

8892608 

5 

2097152 

4408709 

2311557 

60910733 

8569216 

6472064 

70173056 

1 Tetrahedral Meshes 

6 

1337 

2514 

1177 

18729 

4943 

3606 

22642 

7 

10706 

15765 

5059 

81741 

35520 

24814 

164246 

8 

100782 

131204 

30422 

492158 

317367 

216585 

1474817 

9 

220106 

279281 

59175 

956659 

685718 

465612 

3190244 

10 

428538 

533442 

104904 

1694008 

1324614 

896076 

6167300 

11 

2027449 

2452416 

424967 

6818299 

6193783 

4166334 

28862986 


Table 1: Key dehnes the mesh reference; Nb Cells is the number of cells of the mesh; Nb 
dof is the number of discrete unknowns; Nb dof el. is the number of discrete unknowns after 
elimination of cell unknowns; Nb Jac refers to the number of non-zero Jacobian entries after 
elimination of the cell unknowns and equations. 
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HETEROGENEOUS PERMEABILITY - HEXAHEDRAL MESHES 



CPU TIME 


HETEROGENEOUS PERMEABILITY - TETRAHEDRAL MESHES 


IE+00 

lE-01 

lE-02 

lE-03 

lE-04 

lE-05 

lE-03 IE-02 lE-OI IE+00 lE+OI IE+02 lE+03 


VAG Sol 
VAG Grad 



CPU TIME 


Figure 5: Heterogeneous Permeability: Comparison of VAG-FE and HFV on hexahedral and 
tetrahedral meshes. 
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lE-03 lE-02 lE-01 lE+00 lE+01 lE+02 lE+03 

CPU TIME 


ANISOTROPIC PERMEABILITY - TETRAHEDRAL MESHES 



CPU TIME 


Figure 6: Anisotropic Permeability: Comparison of VAG-FE and HFV on hexahedral and 
tetrahedral meshes. 
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Heterogeneous Permeability: VAG 


Hexahedral Meshes 

Key 

Iter 

CPU 

errsoi 

CTT grad 

jump 

^sol 

^grad 

^jump 

1 

8 

1.34E-2 

5.78E-3 

1.74E-2 

8.99E-3 

1.92 

1.97 

1.83 

2 

12 

0.11 

1.53E-3 

4.44E-3 

2.53E-3 

1.92 

1.97 

1.83 

3 

22 

0.98 

3.92E-4 

1.14E-3 

6.72E-4 

1.97 

1.96 

1.91 

4 

41 

8.86 

9.89E-5 

2.91E-4 

1.73E-4 

1.99 

1.97 

1.96 

5 

79 

87.91 

2.48E-5 

7.40E-5 

4.40E-5 

1.99 

1.98 

1.98 

Tetrahedral Meshes 

6 

7 

5.82E-3 

2.01E-2 

0.14 

2.25E-2 

1.80 

0.94 

1.68 

7 

10 

3.73E-2 

5.78E-3 

7.09E-2 

7.03E-3 

1.80 

0.94 

1.68 

8 

20 

0.41 

1.44E-3 

3.52E-2 

1.81E-3 

1.86 

0.94 

1.82 

9 

26 

1.00 

8.11E-4 

2.71E-2 

1.06E-3 

2.20 

1.01 

2.06 

10 

32 

2.11 

5.60E-4 

2.19E-2 

7.36E-4 

1.67 

0.95 

1.62 

11 

53 

12.92 

1.92E-4 

1.31E-2 

2.58E-4 

2.07 

1.00 

2.03 


Heterogeneous Permeability: HFV 


Hexahedral Meshes 

Key 

Iter 

CPU 

errsoi 

CTT grad 

CTT jump 

C^sol 

C^grad 

C^jump 

1 

11 

1.18E-2 

1.34E-2 

4.3E-2 

2.15E-2 

1.94 

1.80 

1.98 

2 

19 

0.13 

3.49E-3 

1.24E-2 

5.44E-3 

1.94 

1.80 

1.98 

3 

35 

1.45 

8.91E-4 

3.41E-3 

1.38E-3 

1.97 

1.86 

1.98 

4 

73 

20.36 

2.25E-4 

9.15E-4 

3.47E-4 

1.99 

1.90 

1.99 

5 

141 

315.38 

5.65E-5 

2.42E-4 

8.69E-5 

1.99 

1.92 

2.00 

Tetrahedral Meshes 

6 

12 

1.56E-2 

l.OlE-2 

0.11 

1.74E-2 

1.88 

0.96 

1.73 

7 

21 

0.22 

2.74E-3 

5.87E-2 

5.24E-3 

1.88 

0.96 

1.73 

8 

43 

3.75 

6.07E-4 

2.75E-2 

1.17E-3 

2.02 

1.02 

2.00 

9 

60 

10.51 

3.38E-4 

2.07E-2 

6.62E-4 

2.25 

1.08 

2.20 

10 

73 

23.52 

2.22E-4 

1.68E-2 

4.37E-4 

1.90 

0.94 

1.87 

11 

119 

166.46 

7.73E-5 

9.87E-3 

1.58E-4 

2.03 

1.02 

1.96 


Table 2: Isotropic test case. Key refers to the mesh dehned in table 1; Iter is the number 
of solver iterations; CPU refers to the solver CPU time in seconds; errsou^rrgrad-, jump are 
the respective L^-errors as defined above; asou ctgrad, ctjump are the orders of convergence of the 
solution, of the gradient and of the jump, respectively. 
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Anisotropic Permeability: VAG 


Hexahedral Meshes 

Key 

Iter 

CPU 

errsoi 

CTT grad 

jump 

^sol 

^grad 

^jump 

1 

7 

6.32E-3 

8.78E-3 

1.98E-2 

8.69E-3 

1.89 

1.99 

1.89 

2 

9 

5.56E-2 

2.37E-3 

4.97E-3 

2.34E-3 

1.89 

1.99 

1.89 

3 

14 

0.67 

6.15E-4 

1.24E-3 

6.06E-4 

1.95 

2.00 

1.95 

4 

26 

6.35 

2.28E-4 

1.57E-4 

3.11E-4 

1.97 

2.00 

1.97 

5 

47 

62.65 

3.95E-5 

7.78E-5 

3.89E-5 

1.99 

2.00 

1.99 

Tetrahedral Meshes 

6 

7 

1.95E-3 

2.73E-2 

0.13 

2.70E-2 

1.95 

0.99 

1.95 

7 

8 

2.14E-2 

7.05E-3 

6.76E-2 

6.98E-3 

1.95 

0.99 

1.95 

8 

15 

0.38 

2.56E-3 

3.92E-2 

2.53E-3 

1.35 

0.73 

1.36 

9 

21 

1.02 

1.34E-3 

2.84E-2 

1.32E-3 

2.49 

1.24 

2.49 

10 

25 

2.24 

9.26E-4 

2.22E-2 

9.14E-4 

1.66 

1.10 

1.67 

11 

41 

13.78 

3.10E-4 

1.36E-2 

3.07E-4 

2.11 

0.95 

2.11 


Anisotropic Permeability: HFV 


Hexahedral Meshes 

Key 

Iter 

CPU 

errsoi 

CTT grad 

CTT jump 

C^sol 

C^grad 

C^jump 

1 

9 

6.02E-3 

2.64E-2 

4.89E-2 

3.35E-2 

1.91 

1.78 

2.01 

2 

16 

8.48E-2 

7.02E-3 

1.43E-2 

8.30E-3 

1.91 

1.78 

2.01 

3 

29 

1.13 

1.81E-3 

3.96E-3 

2.07E-3 

1.95 

1.85 

2.00 

4 

55 

16.55 

4.60E-4 

1.07E-3 

5.19E-4 

1.98 

1.89 

2.00 

5 

108 

248.20 

1.16E-4 

2.86E-4 

1.30E-4 

1.99 

1.91 

2.00 

Tetrahedral Meshes 

6 

10 

1.41E-2 

1.77E-2 

0.14 

1.79E-2 

1.86 

0.98 

1.91 

7 

19 

0.26 

4.86E-3 

7.13E-2 

4.75E-3 

1.86 

0.98 

1.91 

8 

37 

4.56 

1.28E-3 

3.63E-2 

1.21E-3 

1.79 

0.90 

1.83 

9 

47 

12.16 

6.92E-4 

2.62E-2 

6.66E-4 

2.35 

1.25 

2.28 

10 

63 

27.96 

4.75E-4 

2.16E-2 

4.68E-4 

1.69 

0.88 

1.59 

11 

105 

189.66 

1.65E-4 

1.28E-2 

1.58E-4 

2.04 

1.00 

2.09 


Anisotropic Permeability: VAG Lump 


Hexahedral Meshes 

Key 

Iter 

CPU 

errsoi 

CTT grad 

CTT jump 

C^sol 

C^grad 

C^jump 

1 

7 

3.90E-3 

9.09E-3 

2.01E-2 

9.06E-3 

1.89 

1.99 

1.89 

2 

9 

5.15E-2 

2.46E-3 

5.06E-3 

2.44E-3 

1.89 

1.99 

1.89 

3 

15 

0.66 

6.37E-4 

1.27E-3 

6.34E-4 

1.95 

2.00 

1.95 

4 

26 

6.39 

1.62E-4 

3.17E-4 

1.61E-4 

1.97 

2.00 

1.97 

5 

47 

62.19 

4.09E-5 

7.93E-5 

4.07E-5 

1.99 

2.00 

1.99 

Tetrahedral Meshes 

6 

7 

2.11E-3 

2.75E-2 

0.13 

2.73E-2 

1.95 

0.99 

1.94 

7 

8 

2.00E-2 

7.14E-3 

6.76E-2 

7.10E-3 

1.95 

0.99 

1.94 

8 

15 

0.38 

2.60E-3 

3.92E-2 

2.58E-3 

1.35 

0.73 

1.35 

9 

21 

1.02 

1.36E-3 

2.84E-2 

1.35E-3 

2.48 

1.24 

2.49 

10 

25 

2.24 

9.40E-4 

2.22E-2 

9.33E-4 

1.66 

1.10 

1.67 

11 

41 

13.91 

3.15E-4 

1.36E-2 

3.13E-4 

2.11 

0.95 

2.11 


Table 3: Anisotropic test case. Key refers to the mesh dehned in table 1; Iter is the number 
of solver iterations; CPU refers to the solver CPU time in seconds; err soU err grad, err jump are 
the respective L^-errors as dehned above; OisohOigrad,Oijump are the orders of convergence w.r.t. 
of the solution, of the gradient and of the jump, respectively. 
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The test case shows that, on cartesian grids, we obtain, as classically expected, convergence 
of order 2 for both, the solution and it’s gradient. For tetrahedral grids, we obtain convergence 
of order 2 for the solution and convergence of order 1 for it’s gradient. We observe that the 
VAG scheme is more efficient then the HFV scheme and this observation gets more obvious 
with increasing anisotropy. Comparing the precision of the discrete solution (and it’s gradient) 
for VAG and HFV on a given mesh, we see that on hexahedral meshes, the advantage is on the 
side of VAG, whereas on tetrahedral meshes HFV is more precise (but much more expensive). 
On a given mesh, HFV is usually (see [19]) more accurate than VAG both for tetrahedral and 
hexahedral meshes. This is not the case for our test cases on Cartesian meshes maybe due to 
the higher number for VAG than for HFV of d.o.f. at the interfaces F^ on the matrix side. It 
is also important to notice that there is literally no difference between VAG with hnite element 
respectively lumped m/-fluxes concerning accuracy and convergence rate. 


6 Conclusion 

In this work, we extended the framework of gradient schemes (see [7]) to the model problem (4) 
of stationary Darcy flow through fractured porous media and gave numerical analysis results 
for this general framework. 

The model problem (an extension to a network of fractures of a PDF model presented in 
[10], [12] and [3]) takes heterogeneities and anisotropy of the porous medium into account and 
involves a complex network of planar fractures, which might act either as barriers or as drains. 

We also extended the VAG and HFV schemes to our model, where fractures acting as 
barriers force us to allow for pressure jumps across the fracture network. We developed two 
versions of VAG schemes, the conforming hnite element version and the non-conforming control 
volume version, the latter particularly adapted for the treatment of material interfaces (cf. [9]). 
We showed, furthermore, that both versions of VAG schemes, as well as the proposed non- 
conforming HFV schemes, are incorporated by the gradient scheme’s framework. Then, we 
applied the results for gradient schemes on VAG and HFV to obtain convergence, and, in 
particular, convergence of order 1 for ’’piecewise regular” solutions. 

For implementation purposes and in view of the application to multi-phase how, we also 
proposed a uniform Finite Volume formulation for VAG and HFV schemes. The numerical 
experiments on a family of analytical solutions show that the VAG scheme ohers a better 
compromise between accuracy and GPU time than the HFV scheme especially for anisotropic 
problems. 


Acknowledgements: the authors would like to thank TOTAL for its hnancial support and 
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